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C _~i Abstract. In this paper for two-phase parabolic obstacle-like problem, 

An - u t = A+ ■ X{u >0} ~ ^ ■ X{u<0} , (*. x) G (0, T) x Q, 

where T < oo, A+, A > are Lipschitz continuous functions, and O C R™ is 
a bounded domain, we will introduce a certain variational form, which allows 
us to define a notion of viscosity solution. The uniqueness of viscosity solution 
is proved, and numerical nonlinear Gauss-Seidel method is constructed. Al- 
^v^j though the paper is devoted to the parabolic version of the two-phase obstacle- 

like problem, we prove convergence of discretized scheme to a unique viscosity 
I i solution for both two-phase parabolic obstacle-like and standard two-phase 

membrane problem. Numerical simulations are also presented. 

■)—> 

g 1. Introduction 

i— i In this paper we construct and implement a numerical method for the two-phase 

parabolic obstacle-like problem 

^ ( Au - u t = A+ • X{u>o} - A ~ ' X{u<o}, in (0, T) x Q, 

(1) lu(0,x)=g(x), in(0,T)xQ, 
<N [ u (t, x) = h(t), in (0, T) x dn, 

where 

(2) A ± G C 0A (fl) < T < oo, and fl C E" is a given bounded domain. 

' . I Here g{x) is a sign changing continuous function, and the function h(t) is possibly 

^ sign changing function. 

The problem arises as limiting case in the model of temperature control through 
5_i the interior described in |DL761 Section 2.3.2]. 

In the paper SI 'WO!) the authors proved an optimal regularity result for the 
free boundary d{u > 0} U d{u < 0}. They show that if a branch point occurs (i.e., 
the two phases {u > 0} and {u < 0} coexist with vanishing Vw), then nearby 
d{u > 0} and d{u < 0} the boundary is the union of two Lipschitz graphs that are 
continuously diffcrcntiablc in the space variable. 

The stationary case - The two-phase membrane problem has been studied from 
different view points. The optimal C,J C regularity has been proved by Ural'tseva 
[UraOlj in the case when the coefficients A ± are assumed to be constant, and the 
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result was extended by Shahgholian Sha03j for Lipschitz- regular A . The regular- 
ity for the free boundary has been studied by Shahgholian, Ural'tseva and Weiss 
|SW06j . |SUW07| . 

The outline of the paper is the following: In Section 2 we introduce a variational 
form and prove that it has a unique viscosity solution, then using the form we 
construct numerical difference scheme. Convergence of the scheme to the unique 
viscosity solution follows from the so-called Barles-Souganidis theorem. In Section 3 
we develop nonlinear Gauss-Seidel method to solve the discrete two-phase parabolic 
obstacle-like problem, and give some numerical examples of approximate solutions 
by this algorithm. 

2. Convergence of schemes 

2.1. Viscosity solution. We start this section with definition of viscosity solutions 
for parabolic type equations: 

(3) G{D 2 u,Du,u t7 u,t,x) = on Ox(0,T), 

where f2 is a bounded domain and G(M,p x ,p t , u, t, x) is a real- valued discontinuous 
function defined on S™ x R n x R x [0, T] x f2, where S ra is the space oinxn symmetric 
matrices. Here G is always assumed to satisfy the following ellipticity condition 

G(M,p x , Pt ,u,t,x) < G(N,p x ,p u u,t,x) if M > N, 

for all x € R n ,t > 0,Pt e R,p x € R", and M, N G S™ (see [CIL92] ). 

Definition 1. A bounded uniformly continuous function u : [0, T] x ft — > R is called 
a viscosity subsolution (resp. supersolution) for if for all <j> G C 2 ([0, T] x O) and 
all x € Q such that u — <f> has a local maximum (respectively minimum) at (t,x), 
we have 

G{D 2 <j>,D(j>,(j>t,u,t,x) < 0, 

(respectively 

G{D 2 cj),D(j),(j) tl u,t,x) > 0). 

The function u is said to be a viscosity solution of ([3]) if it is both sub- and 
supersolution of For general background about the theory of viscosity solutions 
we refer to |CIL92[ . 

Now we define the following variational equation: 

(4) G[D 2 u, u t , u] = min(u t — Au + A + , max(u t — Au — A~, u)) = 0, 

where A + and A - are positive and Lipschitz continuous as in ([2]). It is easy to see 
that 

G(X,p t ,r) = mm(p t — trace(X) + A + , max(p t — trace(X) — A - , r)) 

satisfies ellipticity and other conditions, as stated in the beginning of this section, 
hence we can apply Definition ([lj as a notion of viscosity sub- and supersolution. 

Lemma 1. If u is the solution (in the weak sense) to (JlJ, then it is a viscosity 
solution to Q. 

Proof. Suppose u* solves the two-phase parabolic equation ([!]) . Then ([!]) will satisfy 
the following inequality in the sense of distributions 

-A" < Au* - u* < A+ in [0, T] x Q, 



A FINITE DIFFERENCE METHOD FOR PARABOLIC TWO-PHASE PROBLEM 



and hence it holds in the viscosity sense as well (see |Ish95j ). We consider three 
different cases: 

• (t,x) G {u* > 0} 

In this case we note that the solution will be a smooth function, in fact 
u* e C x ' t ({u* > 0}), hence one can understand the derivatives in the 
classical sense. Apparently ([!]) will be reduced to 

u* t - Au* + A+ = 0. 

On the other hand 

max(u* — Au* — A~,u*) > 0, 
because of u* (t, x) > 0. Therefore one obtains 

G[D 2 u* 1 u* t , u*} = min(0, max(u* ~ Au* — A~,it*)) = 0. 
Hence u* solves Q as well. 

• (t,x) £ {u* = 0} 

In this case as mentioned above the equation will satisfy the following 
inequality in the viscosity sense 

(5) - \- < Au* -u* t < \ + , in[0,T]xfi. 

Let (t a ,x ) e {u* = 0} and e C 2 ([0,T] x fi), such that u* - 4> has a 
local minimum at (i , ^o)- Then according to ^ and definition of viscosity 
supersolution we obtain 

Acj)(t ,x ) - cf) t (t Q ,x ) < X + (t ,x ). 

Since 

max.(4> t (to,x ) - A(/)(t ,x Q ) - \~(t o ,x o ),0) > 0, 
then we easily obtain 

G[D 2 4>,<t> u u*]{t ,x )>{). 

Thus u* is a viscosity supersolution for our variational equation for all 
(t, x) G {u* = 0}. In the same way if (t , x Q ) e {u* = 0} and ip € C 2 ([0, T] x 
f2), such that u* — ip has a local maximum at (to, Xo), then by taking again 
into account we have 

AV>(£o,a;o) - ipt{to,x ) > -\~(t ,x ), 

and hence we obtain 

G[D 2 iP,iP u u*]{t ,x )<0. 

Thus u* is also a viscosity subsolution at that point. Therefore u* is a 
viscosity solution to Q, for all (t,x) G {u* = 0}. 

• (t,x) G {u* < 0} 

We note here that the solution will be a smooth function on its negativ- 
ity part as well, and u* £ C 2 'l({u* < 0}). Thus one can understand the 
derivatives in the classical sense. 
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In this case ([T]) will be reduced to 

u* - Au* - \- = 0. 
Variational equation Q will lead us to 

G[D 2 u* , u*, u*} = min(u* — Au* + A + , max(0, u*)) 
= min(u* — Am* + A + , 0) 
= min(A~ + A+,0) = 0. 
Hence, in this case u* solves Q as well. 

□ 

Theorem 1 (Uniqueness of viscosity solution). There exists at most one viscosity 
solution of Q. 

Proof. Suppose there exist two viscosity solutions u and v. We have 

(6) mm(u t — Au + A + , max(u t — Au — A~, u)) = 0, 
and 

(7) min(w t — Av + A + , max(w t — Av — A - , v)) = 0, 

in the viscosity sense. 

Denote (It = u(t,x) > v(t, x)}. It is easy to see that 

(8) n T = {u > v > 0} U {u > > v} U {u > > v} U {0 > u > v}. 
since u satisfies we will have for all (j) € C 2 (0 x [0,T]) 

mm((f)t — A(j) + A + , max((/) t — A</> — A~, u)) > 0, 

whenever u — (f> has a local minimum at (t,x). Thus (j) t — Ac/) + A + > 0, therefore 
according to definition ([!]) we conclude that u is a viscosity supersolution to 

u t - Au + A+ = 0, 

at a point (t,x). On the other hand ^ implies also 

mm(ip t — Aip + A + , max(^t — Ai[i — A - , u)) < 0, 

whenever u — ip has a local maximum at (i, a;), and !/j£C 2 (Ox [0,T]). 
By assuming u(t,x) > we conclude 

max(^ f — Aip — A - , u) > 0, 

therefore tp t — Aip + A + < 0. This in turn implies that u is a viscosity subsolution 
to 

u t - Au + A+ = 0. 

Hence we conclude that u is a viscosity solution to u t — Au + A + = 0, for all 
(t,x) e{u> 0}. Thus 

Au — Ut = A + , 
in the viscosity sense whenever u(t, x) > 0. 
Similarly we will obtain that 

Au — Uf = — A~, 

in the viscosity sense, whenever u(t, x) < 0. Again applying the same arguments 
for the cases u > and u < we will get 

-A~ < Au-ut < A+, 
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in the viscosity sense. 

It is apparent that all results will remain true for v as well. Now recalling the 
representation ^ we see that if (t,x) G {u > v > 0}, then by above results 

Au — u t = A + > Aw — vt- 

In a similar way if (t, x) G {u > > v}, we will obtain exactly the same inequality. 
For the cases (t, x) € {u > > v} and (t, x) 6 {0 > u > w} again by above 
arguments we will arrive at the following inequality 

Ait — u t > — A~ = Aw — w t , 

in the viscosity sense. Hence combining all cases together we conclude that 

Au — u t > Aw — Vt, 

for all (t,x) G CI?. Finally we have 

. > f A(w - it) - (u - u) t < in O r , 

[ 1 \ v - u = on <90 T , 

and by applying weak comparison principle for viscosity solutions we obtain 

u(t,x) < v(t,x), for all (t,x) G (It- 

But this inequality is inconsistent with the definition of fix- Hence Ot = 0, and 
for all (t, x) G [0, T] x 11 we have 

u(t, x) < v(t, x), 

and by changing the places of u and v we will obtain the inverse inequality 

v(t, x) < u(t, x), 

for all (t,x) G [0,T] x O. Thus 

7i = V, 

for all (t, x) G [0, T] x il, and this completes the proof of the theorem. □ 

Remark 1. We note that the same approach works for the two-phase membrane 
problem as well. In this case the variational equation will be 

(10) F[D 2 u,u] =min(-Ait + A + ,max(-Au- A", it)) = 0. 

For more information about the two-phase membrane problem we refer to the pa- 
pers |SUW07[ ISW06[ lUraOll IWeiOlj . For numerical analysis we refer to the paper 
[Bozllj .where the author proposed three algorithms. But the numerical methods 
lack of convergence proofs. 

2.2. Convergence for the two-phase membrane problem. Let O be a 

bounded domain in K™, Du and D 2 u denote the gradient and Hessian of u, re- 
spectively and F(x,r,p, X) be a continuous real valued function defined on Q x 
M x K" x S", with 8™ being the space of symmetric n x n matrices. Write 

F[u](x) = F(x,u(x),Du(x),D 2 u(x)). 

Consider the nonlinear, degenerate elliptic partial differential equation with Dirich- 
let boundary conditions, 

l>[tt](aO = infi, 
1 u(x) = g{x) on dfl. 
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We define a uniform structured grid on the domain fi as a directed graph consisting 
of a set of points Xi G CI, i = 1, N, each endowed with a number of neighbors K. 
A grid function is a real valued function defined on the grid, with values Ui = u(xi). 
The typical examples of such grid are 3-point and 5-point stencil discretization for 
the spaces of one dimension and two dimension, respectively. 

A function F h : R N -> R N , which is regarded as a map from grid functions to 
grid functions, is a finite difference scheme if 

F h [u] 1 = F l [u uUl -Ui^-.^Ui- u iK ] (i = l,...,N), 

where {i\, 12, ■ ■ ■ , ik} are the neighbor points of a grid point i. Denote 

F' l [u] = F l [ui,Ui - Uijlj^jjc] = F l [ Ui ,Ui - Uj], i = 1, ...,N, 

where Uj is shorthand for the list of neighbors \ j—Xic- 

Definition 2. The scheme F is degenerate elliptic if each component F l is nonde- 
creasing in each variable, i.e. each component of the scheme F l is a nondecreasing 
function of itj and the differences Ui — Wj . for j = 1, K . 

Since the grid is uniformly structured, we denote h > be the size of the mesh. 
Then in our case the approximation scheme for the elliptic two-phase membrane 
problem is 

(12) F l [m,Ui - uj] = min [Lh,u% + A+, max {L h Ui - A^, m)) , 

where 



^ 1 



(13) L h Ui = 2_^—(ui- u tj ), i = l,...,N. 

3=1 

It is easy to see that F l [iii,Ui — uf\ is non-decreasing with respect to Ui and Ui — 
Uj, therefore the finite difference scheme for two-phase membrane problem is a 
degenerate elliptic scheme. But we know that the degenerate elliptic schemes are 
monotone and stable (see [Obc06 ) . 

Definition 3. (Consistency). We say the scheme F h is consistent with the equation 
(111 at Xq if for every twice continuously differentiable function <p(x) defined in a 
neighborhood of xq, 

F h [4>](x ) -+F[4>](x ) as h^O. 

The global scheme defined on f2 is consistent if the limit above holds uniformly for 
all x € f2. (The domain is assumed to be closed and bounded). 



Lemma 2. The approximation scheme ( 12 ) is consistent. 



Proof. In order to show the consistency we will apply approximation scheme ( 12 ) 
to the twice continuously differentiable function <p(x). 

Suppose xq is a grid point and the function <j>{x) is twice continuously differen- 
tiable around that point. Then if we use Taylor expansion for the function 4>(x) 
around the point xq, we obtain 

K 1 

and 

4>i — > 4>{xq) as h — > 0. 
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Therefore 

max(L h <f>i - \i,<f>i) -> m&x(-A(f>(x ) - X~(x ) ) <f>(x )) > 

and 

+ A+ ->■ -A0(x o ) + A+(x ) 



as ft -> 0. Thus in light of ([10} and ([12} 

F^](x ) -> F[(f>](x ) as ft^O. 

□ 

Now we are ready to formulate the convergence result for the two-phase mem- 
brane problem. 



Theorem 2. (Convergence) The finite difference scheme given by (12) converges 



uniformly on compacts subsets of ft to the unique viscosity solution of the two 



phase-membrane variational equation ( 10 ) 



Proof. By virtue of the so-called Barles-Souganidis Theorem (see |BS91| ) we need 
to show that the scheme is monotone, stable and consistent. The stability and 
monotonicity are provided by Definition |2]), and the finite difference scheme (12 1, 



where we have concluded that it is actually a degenerate elliptic scheme. Consis- 
tency follows from lemma This completes the proof of the theorem. □ 

2.3. Convergence for Parabolic two-phase obstacle-like problem. Define 
Qt = fl x (0,T). As in previous section we consider a uniform structured grid on 
the domain Q consisting of a set of points Xi £ Q, i — 1,...,N, with a number of 
neighbors K. For the time axis discretization we use the following grid tj £ [0, T], 
where j = 1, M. 

Unfortunately the notion of degenerate elliptic schemes is not applicable to this 
case, since the variational form Q is not degenerate parabolic. To be more clear 
it is noteworthy that the degenerate elliptic schemes defined above and applied for 
elliptic version of the Two-Phase Obstacle-like problem, are just particular case of 
the schemes defined in the Barles-Souganidis result. But since the variation form 
of Parabolic two-phase obstacle-like equation is not degenerate parabolic as stated 
above, then we have to proceed all the steps to check whether our scheme satisfies 
the required conditions for monotone schemes stated in the Barles-Souganidis the- 
orem or not. In order to do that we follow the notations of [BDR95 . A numerical 
scheme can be written as 

S(m, u) = S(At, Ax, m, j, it™, u) = 0, 

for 1 < j < N and 1 < m < M, where N and M are respectively the number of 
grid points in space and in time. Here u denotes the vector (ui)fc / in IR^* 1 . Finally 
At and Ax denote the time and the space mesh size respectively. 
The definition of Monotonicity for the scheme will be as follows: 

S{At,Ax,m+l,j,uJ +l ,u) < S{At,Ax,m+l,j,vJ- +1 ,v) \fu>v and if uf +1 = v™ +1 , 

for any At, Ax > 0, 1 < j < N, 1 < m < M, and for all u and v in R NM . In our 
case the scheme will have the following form 

(14) S(m + l,u) = min(5(u) + AtX+j; max^u) - AtA - j, Atuf + 1 )) = 0, 
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where 



At - 
S(u) = u™ +1 - uf + t^Lu™, and Lu? = ^ « - u£), i = l,...,N, 

^ ' q—l 



here {ii, . . . , ix} are the neighbor points of a grid point 



Lemma 3. The scheme (14) is monotone and stable provided the following non- 
linear CFL condition holds 

< 15 > £ 



Proof. In order to prove monotonicity of ( 14 ) it is enough to see that 

S(u) < S(v), iiu>v and if uf +1 = vf +1 , 

for any At, Ax > 0, 1 < j < N, 1 < m < M, and for all u and v in R NM . 
We have 

§(u) = u T ^u T+ ^ ek-o = < l -«r (i - ^) - ^ g 

Since w' n+1 = and > vf for all 1 < j < N, 1 < m < M and from CFL 

condition 

A< 

(Axf ~ 1 ' 

we arrive at 

S(u) < S(v). 

To obtain the stability we refer to the Lemma 4.1 in |UL.I09 |, where the authors 
proved comparison principle of numerical scheme defined for one-phase parabolic 
type equation arising in American option valuation. It is not hard to see that we 
can do the same induction for our scheme as well (this is standard). Once we have 
this, the stability follows directly, because we can mimic with the boundary values 
of the scheme and see that for fixed spatial and time discretization our scheme will 
stay between the maximum and minimum values of the discrete boundary values 
which are fixed a priori. This part is also standard to proceed for such schemes, 
that's why we skip the detailed proof. 

Thus we have the stability and monotonicity for ( 14 ) . □ 



Consistency of (14) can be done as in previous section for the two-phase mem- 
brane problem. We can easily observe that the following limit holds 

5(Ai,Ax,m + l,j,<?!>™ +1 >) 
(16) At J ' G[D^, 0, 

as At, Ax — > 0, for every twice continuously differentiable function <j){t,x). 

Remark 2. It is easy to see that 3— point and 5— point stencil discretization will 
lead us to the following CFL conditions 

At 1 , At 1 

< - and j-r — r-^- < — , respectively. 



(Ax) 2 " 2 (Ax) 2 ~ 4' 

Remark 3. We note that we can consider the implicit discretization of this scheme 
as well. In this case we will have unconditionally monotone and stable scheme. 
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We are ready to write down the main result of this section. 



Theorem 3. (Convergence for parabolic case) The solution u o/( 14 ) converges as 
At, Ax —> uniformly on compacts subsets of Qt to the unique viscosity solution 
of the two-phase parabolic obstacle-like variation equation Q . 

Proof. The proof is again immediate consequence of Barles-Souganidis theorem as 
in previous section. The stability and monotonicity are provided by lemma ([3]). 



The consistency follows from the limit (16). □ 



3. Numerical method and Simulations 

3.1. Numerical method. As mentioned in the abstract of the paper for construct- 
ing a numerical method we refer to Nonlinear Gauss-Seidel method. 
Suppose u m is a shorthand of (iff")j. We proceed as follows: 

• First Step. 

u m +i = min (u m - -^-r Lu m + At A - , 

\ (AxY 

• Second Step. 

u m+1 = max ( u m - -^-Lu m - AtA+, u m ' 
(AxY 



In order to see the consistency of the method with the difference scheme ( 14 ) , we 
eliminate u m+ 2 in the above equality. We obtain 

u m+1 = max (u m - -J^-Lu™ - AtA+,min ( u m - -^-Lu m + AtA~,0 
\ (Ax) 2 \ (AxY 

Therefore 

min ( u m+1 -u m + -£L-Lu m + At\ + , max ( u m+1 - u m + -^L-Lu m - AtA", ) ) = 0. 

\ (AxY V (Aa;) 2 

Dividing the first argument in the max and in the min by At we will derive the 
desired consistency condition. 

3.2. Simulations. In this section we present some numerical examples for the two- 
phase parabolic obstacle-like problem. For all examples we consider £1 = [0, 1] and 
T =1. 

Our equation reads: 

Uxx -ut = X+X{u>o} - A~X{k<o}, (t, x) e (0, 1) x (0, 1), 
(17) { u(0,x)=g(x), 

u(t,0) = hx(t),u(t,l) = h 2 (t). 



In the figures T][3 for different A + and A numerical simulations are shown. For 



all cases we take the initial data at time t = to be linear and the boundary values 
h\(t) and ft.2(t) 1° be constant. We use implicit discretization in space and forward 
discretization in time. Numerical examples were constructed with 200 discretization 
points in space and 250 discretization points in time. In the contour plots of figures 
are clearly visible the positivity and negativity sets of solutions. In all cases we 
observe that as time evolves, the free boundary becomes more stable and after 
some amount of time it does not change much. This is expected since from the 
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theory of Parabolic Two-phase Obstacle-Like problems we know that ||wt||oo — > 0, 
as t — > oo. 





4. Conclusion 

In the paper we apply the theory of viscosity solutions to the Parabolic Two- 
Phase Obstacle-Like problem in order to develop a convergent numerical scheme. 
We have observed that the method is applicable for the elliptic case as well. In 
section 3 we applied the numerical method for different X ± , and concluded that 
the developed theory works well. One of the main questions remain to consider the 
rate of convergence of the above scheme, which is the subject of the forthcoming 
paper by the author and his collaborators. 
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(a) Numerical solution (b) Contours 

FIGURE 3. The left picture shows a numerical solution for A + = A~ = 
0.6 and g(x) — 8x — 4. The right picture shows the contours of u(t,x). 
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